Numerical study of vortex matter using the Bose model: 
First-order melting and entanglement 
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We present an extensive numerical study of vortex matter using the mapping to 2D bosons and 
Path-Integral Monte Carlo simulations. We find a first-order vortex lattice melting transition into an 
entangled vortex liquid. The jumps in entropy and density are consistent with experimental results 
on YBa2Cu307_,5. The liquid is denser than the lattice and has a correlation length l z ~ 1.7eao 
in the direction parallel to the field. In the language of bosons we find a sharp quantum phase 
transition from a Wigner crystal to a superfluid, even in the case of logarithmic interaction. We also 
measure the excitation spectrum of the Bose system and find the roton minimum to be insensitive 
to the range of the interaction. 
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I. INTRODUCTION 

The Abrikosov mean-field theory of the vortex system 
predicts a peculiar continuous freezing transition, involv- 
ing two symmetries, at the upper critical field H c #i. The 
gauge symmetry is broken in the direction parallel to 
the field by the appearance of a superconducting order 
parameter \&(x). The presence of the vortex lattice, on 
the other hand, leads to a periodic modulation of the 
amplitude of the order parameter, breaking the transla- 
tional invariance of the system. This behavior contra- 
dicts a general argument by Landau, stating that break- 
ing of continuous translation symmetry should be asso- 
ciated with a first order transitionu. In order to resolve 
the paradox, we have to go beyond mean-field theory 
and take thermal fluctuations into account: The fluctu- 
ations melt the vortex lattice in a first-order transition 
at temperatures below the mean-Jjield transition in agree- 
ment with the Landau argumcntlj. The size of the region 
where thermal fluctuations are important is determined 
by the Ginzburg criterion: In conventional superconduc- 
tors this region is small and mean-field theory offers a 
good description of the relevant physics. In the high-T c 
materials, on the other hand, thermal fluctuations drasti- 
cally change the mean-field phase diagram. Understand- 
ing the effect of thermal fluctuations is therefore vital for 
the understanding of the vortex state in high-T c super- 
conductors. 

There exist a number of fundamental questions asso- 
ciated with the melting of the vortex lattice. To begin 
with, the existence of two symmetries naturally leads to 
the question of whether these have to be broken simulta- 
neously in one transitionQ. In terms of vortices, the gauge 
symmetry along the field is destroyed by the alignment of 
the vortices with the applied magnetic field, whereas the 
translational symmetry is broken by the ordering of the 
lines in the direction transverse to the field. A melting 
scenario involving a disentangled liquid, which has been 



the tapic of an intense scientific discussion over the past 
ycaraj'u, would break these symmetries in two separate 
transitions. A further topic deals with the properties of 
the melted phase: Is the resulting vortex liquid thermo- 
dynamically identical to the normal state or do further 
phase transitions exist at higher temperatures. 

Although the importance of thermaLfkictuations in the 
vortex system was recognized earlyLrBa, most progress 
has been made only in recent years. One important 
reason for this is that it has now become possible to 
observe the vortex lattice melting transition experimen- 
tally: After the first indirect observations based on resis- 
tivity measurementserL2l, more recent experiments have 
observed the first-order transition directly by .measuring 
a jump in the magnetization or the latent heatEJIi-l. The 
very large latent heat which was measured in many ex- 
periments was puzzling at fcsfc^but the problem has now 
been resolved satisfactorilytSta 

On the theoretical side, there still does not exist a 
reliable theory for the vortex lattice melting, transition. 
Early work usingr-the renormalization grouptl or density 
functional theoryO have indicated a first-order transi- 
tion. Elastic theory has been combined with the Lindc- 
mann criterion to produce a melting .line, in good agree- 
ment with experimental observationsEallj. Even though 
the Lindemann criterion does not give explicit informa- 
tion on the character of the melting line, it can be argued 
that the fact that it works so well gives support for a 
first-order transition. The question regarding the possi- 
ble existence of a disentangled liquid has been addressed 
by analyzing the energy and stability of entangled con- 
figurations both in the liquid and in the lattic<£3, with 
the conclusion that no disentangled liquid should exist. 

The absence of an analytical description of the vortex 
lattice melting transition has led to a large interest in 
numerical simulations- A-,popular approach is to use the 
frustrated XYjaod-elEirEl or the closely related Lattice 
London mode&B. Unfortunately, many of the simu- 
lations have suffered from highly non-trivial finite-size 
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effects and have indicated two transitions instead of one 
at low filling factors, whereas a first-order transition was 
seen at large filling, factorsE3. The problems have been 
overcome recentlycJ, and a picture with a single first- 
order transition is emerging. An overview of these prob- 
lems can be found in Ref. |26|. Simulations based on the 
Lowest Landau Level (LLL) approximation have shown 
a first-order vortex lattice melting transition^ and pb- 
tained entropy jumps in agreement with experiments!!!! 

In this work we adopt a different approach from those 
described above. It was suggested by Nelson that the vor- 
tex system is equivalent to a system of interacting bosons 
in two dimensions^. This idea allows for using standard 
many-body techniques in describing the vortex system 
and has bqc^pif a popular approach for analytical work 



on vorticesE3tII. The mapping to bosons, which we shall 
refer to as the Bose model in the following, has not been 
widely used for numerical work, however. This is perhaps 
surprising, since there exist numerical algorithms for sim- 
ulating interacting bosons exactly. The main advantage 
of the Bose model lies in its generality: By analyzing 
this model, which contains very few adjustable param- 
eters, we are able to obtain results for vortices, bosons, 
and interacting elastic lines, such as polymers, in gen- 
eral. Our approach also allows us to answer two questions 
which cannot be considered with the methods above. The 
first is obvious: How well does the Bose model describe 
the real vortex lattice? This has not been studied be- 
fore, as it is difficult to make any quantitative predictions 
for strongly interacting bosons analytically. The second 
question is of a more fundamental nature: It has been 
known for a while that the 2D Bose Coulomb Liquid-d 
not have a Bose-Einstein condensate even at T = 0E3~ 
Consequently, it has been questioned whether this sys- 
tem necessarily has to be crystalline or superfluid in i, 
ground state, or whether some other phase might existl 
Here we show that the system has a direct transition from 
a crystal to a superfluid, and that a superfluid without a 
Bose-Einstein condensate can indeed exist. 

Some early results of this work have been published 
previouslytid. Here we expand on our previous work, 
adding new results for melting in a compressible system 
using an isobaric Monte Carlo method. These results 
allow us to use the Clausius-Clapeyron relation as an 
consistency check on the simulation, since the change in 
density and the change in entropy can be measured inde- 
pendently at the transition. New results for the energy of 
corresponding Bose system are presented, and we discuss 
all results both in the language of bosons and vortices, 
further emphasizing similarities and differences between 
the two systems. Finally, specific details on the simula- 
tion are provided. 

The outline of the paper is as follows: In Sec. || we 
describe the model and discuss its applicability to the 
vortex problem and the approximations made. We then 
turn to a qualitative description of the bosonic phase- 
diagram in Sec. [II, and relate it to the phase diagram for 
vortices. The algorithm and the simulations are discussed 



in Sec. IV. Numerical results for the 2D Bose system are 
given in Sec. ^ and the numerical results for the vortex 
lattice melting transition are presented in Sec. [v|. In 
Sec. ^ 
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we look at some properties of the vortex liquid 
and our conclusions are summarized in Sec. VIH| . Three 
appendices provide technical details on various aspects 
of the simulation. 



II. DESCRIPTION OF THE MODEL 

A. Boson-vortex mapping 

We consider a system of N bosons in two dimensions 
interacting through the potential V(r) — g 2 Ko(r/A), 
where g 2 is the energy scale of the interaction, Kq(x) 
is a modified Bessel function of the second kind, and A 
measures the range of, the interaction. In the Feynman 
path-integral pictural3, the system is described by in- 
teracting "elastic" world-lines with the imaginary-time 
action 
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Here, the two-dimensional vectors Rj(r) represent the 
positions of the bosons as a function of the imaginary 
time r and T is the temperature of the system. The 
boundary condition for bosons is Ri(0) = Rj(/3), i.e., 
every line ends either on itself or on some other line. 
The partition function is the average of all possible line 
configurations subject to this boundary condition and 
weighted by the action (0): 



Z(T,A,N) ee e-nwyT = -1 /"fjDR, 
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where F(T, A, N) denotes the Helmholtz free energy and 
A is the area of the system (we are considering two dimen- 
sions) . It is often useful to consider a system of bosons 
at fixed chemical potential /i rather than at fixed number 
of particles. Switching ensemble is easily done according 
to 

OO 

Z(T, A, (j,) ee e -n(T,A,»)/T = Z (T, A, N)e^ T , (3) 

in which case f2(T, A, /x) denotes the grand-canonical po- 
tential. 

It was pointed out by Nelson, that the action in 
Eq. (0) can also be interpreted as the London free- 
energy for a system of interacting vortices in type-II 
superconductorsOLJ. The starting point is the London 
free-energy functional 
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where <£>o is the flux quantum, £o = ^oA^ 71 "-^) 2 is the 
vortex line energy, L z is the thickness of the sample, and 
ei is the elasticity of a vortex line. The latter depends 
strongly on the anisotropy and we have the estimate 
ei ~ £ 2 £o ln^o^y 7 ^), where e 2 < 1 is the anisotropy 
parameter, £ is the coherence length and do is the lattice 
spacing . A detailed discussion of the line energy is given 



in Sec. [IB 



We will refer to Eq. fl|) as the Bose model 
for the vortex system in the following. 

Comparing Eqs. (0) and (||), and requiring T/T = 
S/h, we obtain the following formal mapping between 
2D bosons and vortices: 
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H/Tb <-> L z 
m <-> Ei 



(5) 



The temperature of the Bose system is denoted by T B , 
in order to distinguish it from the (physically different) 
temperature T of the vortex system. Some conclusions 
can immediately be drawn from this mapping: The tem- 
perature of the vortex system corresponds to the Planck 
constant of bosons, implying that thermal fluctuations of 
vortices are equivalent to quantum fluctuations of bosons. 
Quantum fluctuations are more relevant for bosons with 
small mass. Since the mass of the bosons corresponds 
to the elasticity of the vortices, it follows that thermal 
fluctuations are more important for vortices with small 
elasticity, i.e., for anisotropic materials. Finally, the vor- 
tex system in a bulk superconductor, with L z — ► oo, cor- 
responds to the ground state, Tb = 0, of the bosons. 

As in the case with bosons, it is often useful to consider 
a system where the number of vortices is allowed to fluc- 
tuate. The number of vortices determines the magnetic 
field and we therefore have to consider the Gibbs energy 
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This quantity should be compared with S/h— Nfi/Ts for 
the bosons, and we obtain the formal equivalence 
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The chemical potential fi is negative for H < H c i — 
$o in k/47tA 2 , in agreement with the fact that no vortices 
enter the sample for fields below the lower critical field. 

It will be convenient to introduce dimensionless param- 
eters in the following analysis, which is done by choosing 
a length and an energy scale. As the length scale we take 
the lattice constant ao of a triangular lattice, i.e., the 



density is p — 2/<ZqV3- The natural choice for the energy 
scale is g 2 for bosons and Eo^o for the vortex system. If 
we consider a fixed number of particles, where the chem- 
ical potential is not needed, the thermodynamic proper- 
ties of the system are completely determined by the three 
dimensionless parameters 
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with the dimensionless action 
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(9) 



The parameter A is the de Boer parameter and measures 
the size of fluctuations; quantum for bosons and thermal 
for vortices. This is the most important parameter in the 
present study. (3 is the dimensionless inverse temperature 
and we will only be considering the limit j3 — * oo here, 
corresponding to the ground state for bosons or vortices 
in a bulk superconductor. Finally, it is often a good ap- 
proximation to use A = oo, which further reduces the 
number of parameters. In the latter case, we are study- 
ing the two-dimensional Bose Coulomb liquid (2DBCL) 
with logarithmic interaction. Note that simply taking the 
limit A — * oo in the action (Q) leads to a diverging energy 
and the interaction with a background charge therefore 
has to be subtracted. The density is then fixed by the re- 
quirement of charge neutrality and the system is therefore 
incompressible. In order to have a finite compressibility, 
which is needed for a density (or magnetization) jump 
at a first-order phase transition, we have to use a finite 
value of A. Since the vortices in the XY-model interact 
logarithmically, it follows that this model will neveiuhave 
a magnetization jump, in contrast to recent claimso. 



B. Validity of the model 

Before discussing the phase diagram of the model de- 
scribed above, we want to discuss its validity as a model 
for vortices in a type II superconductor. It is natural to 
divide this discussion into three parts, considering bound- 
ary conditions, linearization, and retarded interactions 
separately. 

The Bose model differs from a real vortex system in 
the choice of boundary conditions: The natural boundary 
conditions for vortices are given by 
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since the screening currents, which encircle the vortex, 
cannot leave the sample. Since one cannot hope to simu- 
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late a sample with realistic dimensions we make use of pe- 
riodic boundary conditions in all directions, as this min- 
imizes the effect of the boundaries by ensuring V • B = 
everywhere. The expectation then is that the extrap- 
olation of our numerical results to infinite system size 
provides an accurate description of the transition. The 
bosonic boundary conditions in the longitudinal direc- 
tion, Ri(/3) = Rj(0), correspond to periodic boundary 
conditions for the magnetic field in the z-direction and 
have been used in most numerical simulations of vortices 
in bulk materials so far. 

Another question is related to the linearization leading 
to the elastic term in Eq. (0), i.e., 
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for an isotropic superconductor. This expansion is only 
valid as long as the vortex lines do not deviate too much 
from the z-axis. It is very accurate, however: For a vor- 
tex at a 30° angle with the z-axis, the error is roughly 1%. 
Close to the melting transition we have (dH/dz) 2 « c? L , 
showing that the Taylor expansion is valid. Scaling the- 
ory proves this, result to be true also in an anisotropic 
superconductoroEl 

The use of the z-coordinate as a parameter for the 
vortex lines prohibits the formation of vortex loops in 
the a&-planes. Since these loops have-been argued to be 
important for vortex lattice meltingrHrl, we would like 
to remark on this point here. The free energy for a planar 
vortex loop of length L is given by 



F L 



(12) 



if we use a simple random walk argument for estimating 
the entropy. Loops will begin to proliferate when the free 
energy becomes negative, i.e., for T > ££o£/ m 3. This, 
however, is essentially the Ginzburg criterion, showing 
that planar vortex loops are important in the fluctuation 
regime close to the upper critical field H C 2- Due to the 
smallness of the Lindemann number for vortex lattice 
melting, the melting Una. is usually located well outside 
the fluctuation regimdj'Ej. 

A more important approximation is the way the Bose 
model, i.e., Bosons with Yukawa interaction, neglects the 
"retarded" interaction of the real vortex lattice by the 
substitutionQ 



-^R l3 (z.z>y+(z-z>yi\ 
dz' , . = . = ^K [R ij (z,z)/X}. (13) 



^R l3 {z,zf + {z-z'f 



We have used the notation i?y(z, z') = |R^(z) — R.,-(z')|. 
Again, this approximation is valid as long as the lines 
are straight, at least over distances A in the z-direction. 
However, Eq. (|l3|) replaces a non-local interaction in the 
z-direction with a local one. The consequences of this 
can be understood in terms of the elastic properties of 



the vortex lattice: Since it is the non-local character of 
the interaction which produces dispersion, it follows that 
the substitution ( |l3| ) will remove any fc z -dependence in 
the elastic moduli. The shear modulus cqq is not disper- 
sive and remains unchanged. The compression modulus 
c\\, which is not important for vortex lattice melting, is 
changed in a trivial way. The tilt modulus, however, de- 
serves a more detailed treatment. It is easy to see that 
the tilt modulus corresponding to Eq. M) is given by 



c 44 



(14) 



because only the self-energy terms contribute to the 
tilt energy. The real tilt modulus of the vortex lattice 
is rather |-e<j«Hpljcated, especially in the case of layered 
materialscaOES. The single vortex contribution is given 
by 
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(15) 



and has the same form as Eq. ( |l4| ) apart from a logarith- 
mically weak fc z -dispersion. In addition, the real vortex 
system has a bulk contribution to the tilt, emanating 
from the interaction between vortices, which is given by 



cu(K, k z 



B 2 



4tt l + X 2 K 2 /e 2 + X 2 k 2 



(16) 



This term is absent in the Bose model, since the interac- 
tion energy only depends on the distance between lines 
and not on their angle with the applied field. For ther- 
mal fluctuations, the relevant modes have wave vectors 
K K BZ k, ^/An/ao and k z ~ K/2e. Whence, the 
single- vortex part of the tilt modulus gives the dominant 
contribution to the stiffness of the lattice and we expect 
the Bose model to be valid for describing thermal fluctu- 
ations as long as we use Eqs. (ful) and (15J) to determine 
£l, i.e., by choosing si = e 2 Eo \ma /2y / Tr^). 

The above discussion made use of the elastic moduli of 
the vortex lattice and is not directly applicable to the liq- 
uid. We do not know enough about the properties of the 
vortex liquid to be able to make an analytical comparison 
between the Bose model and the London model in this 
case. A major difference is that the vortex lines in the 
Bose model always repel each other, whereas the inter- 
action between real vortices is proportional to the cosine 
of the angle between the lines and can become attrac- 
tive. The difference is not important in the lattice, since 
the vortices are almost parallel anyway. In the liquid the 
vortices are not so well aligned and the difference will be 
larger. We argued above that the linearization of the self 
energy is very accurate, since the error is a fourth-order 
correction in dRi/dz. Similarly, the error done by drop- 
ping the cosine is a second-order correction in the same 
quantity. The Bose model therefore still contains the 
main part of the interaction between the vortices and we 
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expect our results to be in rough quantitative agreement 
with those of the real system. This is further supported 
by the fact that our model correctly captures the prop- 
erties of the vortex lattice melting transition, as will be 
shown below. 

A very interesting consequence of the direction- 
dependent interaction for real vortices is the existence of 
an attractive van der Waals interaction even for straight 
vortices. This attractive interaction has important con- 
sequences for the low-held phase diagram of anisotropic 
superconductors, below the lower melting lindlj. For 
the magnetic inductions considered here, however, with 
ao < A, the van der Waals interaction can safely be ig- 
nored. 



III. THE PHASE DIAGRAM FOR BOSONS 
WITH LONG-RANGE INTERACTION 

The phase diagram for Coulomb interacting bosons 
comprises three phases and is shown schematically in 
Fig. |l|. The parameter A measures the size of quan- 
tum effects, and the line A = therefore corresponds 
to the classical 2D Coulomb plasma. As the temper- 
ature is lowered, (3 is increased and we expect a tran- 
sition from a fluid to a crystal. Numerical simulations 
agree on tius-transition being hrst order, taking place at 
(3 m w 140cJr2l. For large quantum effects, i.e., for large 
A, the low temperature phase is a superfluid. The tran- 
sition from a classical fluid to a superfluid is generally 
accepted to-be a Kosterlitz-Thouless transition in two 
dimensionsEd Finally, at large /3, the system goes from 
a lattice to a superfluid as the size of quantum fluctua- 
tions is increased. Since no temperature is involved, this 
is a quantum phase transition- and the critical value of 
A is A m 1=3 0.062 for /3 — oolj. This transition, which 
is equivalent to the vortex lattice melting transition in a 
superconductor, is the topic of the present study. 

The phase diagram can be conveniently discussed in 
terms of three energy scales (potential, kinetic, and ther- 
mal) and two symmetries (gauge and translational). As 
long as the thermal energy Tb = g 2 / dominates, the 
system possesses both global gauge invariance and con- 
tinuous translational symmetry. If the potential energy 
g 2 dominates at low temperature, the translational sym- 
metry is broken and we find a crystal. Otherwise, if the 
quantum mechanical kinetic energy g 2 A 2 dominates, the 
gauge symmetry is broken and the system becomes super- 
fluid. Strictly speaking, since we are in two dimensions, 
none of the symmetries is really broken at finite tempera- 
tures. Rather, both the crystal and the superfluid phase 
are characterized by quasi long-range order, as will be 
discussed below. 

According to the phase diagram in Fig. [j] two sym- 
metries must change simultaneously at the phase bound- 
ary between the crystal and the superfluid. This has led 
to speculations on the existence of other phases, such 



as the supersolid, in which both symmetries are broken. 
The first suggestion for the existence of such a super- 
solid phase Sn the context of superfluid He II was made 
by AndreevE], see Ref. for an overview. The topic has 
also been studied in the context of Josephson junction ar- 
rays, where a supersolid has been observed numerically^], 
and for voxtices in superconductors, where no supersolid 
was foundElj. 

Another suggestion is due to the rather peculiar prop- 
erties of the 2D Bose Coulomb Liquid. It can be shown 
that the off-diagonal long-range order in the ground state 
decays as 

(^(r)#))"-°, (17) 
where the exponent a is given 

l/VsV 72 / m 5 2 \ 1/2 ^0.186 
a= A[l^) =U^J *— • (18) 

The system therefore does not have true off-diagonal 
long-range order (ODLRO) even in the ground state. 
This is in contrast to the case with short-range inter- 
actions where ODLRO decays algebraically only at fi- 
nite temperatures-according to the Hohenberg-Mermin- 
Wagner theoremEll. The effect can be understood in 
terms of the long-range interaction suppressing long 
wavelength density fluctuations. Because of the uncer- 
tainty principle, this implies that the fluctuations in the 
conjugate variable, the phase, have to diverge. 

The question is how large a can become before super- 
fluidity is destroyed. One estimate is based on the anal- 
ysis of the momentum distribution function: Computing 
the Fourier transform of Eq. ( [l7j ) we obtain 

n q oc q a ~ 2 , (19) 

and the momentum distribution function actually van- 
ishes for q — > if a > 2. Since the crystal melts at a m « 3 
(A m pa 0.062), it is relevant to ask whether a 2D Bose liq- 
uid with 2 < a < a m can be superfluid. Otherwise, a 
non-superfluid Bose liquid would exist at Tb = 0. One of 
the results of this work is to show that a can indeed be 
larger than 2 in a superfluid, and that superfluidity only 
disappears when the lattice forms, confirming the conjec- 
ture that a Bose system is either a superfluid or a crystal 
at Tb = 0. Note that the concept of a superfluid with- 
out a Bose-Einstein condensate is really not as striking 
as it might seem: A 2D Bose system does not have a real 
condensate at finite temperatures but can still be super- 
fluid. Also, for He II in three dimensions, the condensate 
fraction is below 10% at zero temperature- whereas the 
superfluid density equals the total density^. 

A further consequence of the long-range interaction is 
that the solid phase appears at low densities, i. e., it is 
a Wigner crystal. For a finite interaction range A, we 
also find a superfluid phase at low densities with A < 
oq. Finally, we note that since the interaction is purely 
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repulsive no generic liquid exists; the state which we refer 
to as a liquid might equally well be called a gas or a 
plasma. 

It is instructive to discuss the phase diagram directly in 
terms of world lines. The partition function (^) includes 
all particle permutations, i.e., all ways of connecting the 
lower ends of the world lines at t = with the upper ends 
at r = P, and most of the possible line configurations are 
heavily entangled. There are two factors that prevent 
this entanglement: The length of the world lines is de- 
termined by the temperature of the Bose system and the 
world lines can therefore only entangle if the temperature 
is low enough. A simple estimate for the entanglement 
temperature is given by 
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(20) 



marking the appearance of superfluidity. Since the in- 
teraction between the lines is also less effective at high 
temperatures, it follows that the high-temperature phase 
is a normal liquid. If the interaction is strong enough, 
it will prevent entanglement even of infinitely long world 
lines, and we have a crystal ground state. 

The above phase diagram can be reinterpreted as a 
phase diagram for vortices by use of the parameter map- 
ping in Eq. (^). The Wigner crystal corresponds to the 
vortex lattice, the supcrfiuid is equivalent to an entan- 
gled vortex liquid, and the normal fluid corresponds to 
a disentangled vortex liquid. Since the parameter /3 is 
proportional to the thickness of the sample, we are only 
interested in the limit (3 — > oo for a bulk superconductor. 
It is important to use a large enough L z in a numerical 
simulation: The dashed line in Fig. ^ shows the curve 
obtained by keeping B (ao) constant and varying T for 
some finite L z . The curve passes from a lattice, via a 
disentangled vortex liquid, into an entangled liquid. As 
L z is increased, however, the curve is pushed to higher 
values of /3 and we have a direct transition from a lattice 
to an entangled liquid. 



IV. NOTES ON THE ALGORITHM 

One of the advantages of using the Bose model for 
describing the vortex system is that there exist well es- 
tablished algorithms for simulating bosons numerically 
exactly. The ground state of a Bose system can be stud- 
ied using variational and diffusion Monte Carlo methods, 
which has been done for, the system considered here by 
Ceperley and MagrcE3o. They observed the disappear- 
ance of the Wigner lattice as A was increased and also 
studied the momentum distribution function in the liq- 
uid, which agreed with the result in Eqs. (17) and (18). 
Simulations at Tg = are usually not very good at com- 
puting response functions, such as the superfluid density, 
since these depend on the excitations available in the 
system rather than on the ground state wavefunction. 



Although this difficulty can be overcome for the special 
case of the superfluid densityEj, we have chosen the more 
conventional approach of simulating a system at finite 
temperature. Not only does this allow us to compute 
the superfluid density using the winding number, but it 
also allows us to study the excitation spectrum from the 
dynamic structure factor. 

Exact finite temperature simulations on Bose systems 
were introduced by Ceperley and Pollock and aualied 
to superfluid He II in two and three dimensionsEjEj. 
The method is known as the Path Integral Monte Carlo 
(PIMC) method, and differs from the more widely used 
world-line methods in that it does not introduce a lattice. 
It is therefore possible to obtain numerically exact results 
for realistic interactions between the bosons, as has been 
amply demonstrated by Ceperley and co-workers in com- 
puting the superflui|d-idensity and the excitation spectrum 
for superfluid He HE3. 

The method involves discretizing the imaginary time 
in the path integral for the partition function, which, in 
the Trotter approximation, becomes 



Z((3,A,N) 



with 



1/1 



N\ \2itA 2 t 



MN „ M N 



nn^ 



m— 1 i—1 



S[{Ri]} 



(21) 



■S-{R;H=y: (R<,m+ . 1 A ,. Ri|m) + E -Ko(iWA). 



i.m 



2A 2 7 



(22) 



Here the indices label the particles, m labels the 

Trotter slices ranging from 1 to M, and we have r = j3/M 
for the size of the imaginary time step. We have typically 
used M = 100 and N = 25 — 100 in the simulations pre- 
sented here, with periodic boundary conditions in all di- 
rections. A rhombically shaped system, with the smaller 
angle equal to 60°, has been used to allow for a trian- 
gular lattice without frustration. Since the range of the 
interaction is usually larger than the size of the system, it 
is important to account for the periodicity of the system 
also in the interaction. We do this by solving the London 
equation with the correct boundary conditions and the 
resulting expressions are given in Appendix |A|. 

Computation of thermodynamic averages using Eq. 
( pl| ) involves evaluating a multi-dimensional integral nu- 
merically. Since the dimensionality is high, with MN ss 
10 3 — 10 4 , an efficient algorithm is needed. Most of the 
PIMC algorithm, including the bisection method for gen- 
erating new-Line configurations, is described in the review 
by CeperleyEj. There is one point, though, where our al- 
gorithm differs from those used before: In order to cap- 
ture the quantum effects correctly in the Bose system, 
it is necessary to sample the particle permutations ac- 
curately. Thus we have to probe all ways of connecting 



G 



the lower ends of the lines Ri(0) with the upper ends 
Rj(/3). Usually this is done by allowing a small number 
of lines to change their endpoints in every Monte Carlo 
update. We were unable to equilibrate the system using 
any such algorithm. Rather we use a random walk algo- 
rithm which allows the system to choose a permutation of 
a large number of lines according to its statistical weight. 
The algorithm is described in Appendix ||. For N = 49 
or less, we actually do this random walk in the space of 
the permutations of all N lines, leading to a reasonably 
fast convergence of the winding number (see below). 



A. Thermodynamic properties 

We now turn to the measurable quantities of the sys- 
tem. The existence of a lattice is tested by measuring 
the structure factor 5(Q), 

1 



S(Q) = -<p(QM-Q)> 



MN \E 



cxp {iQ ■ (R 4 , m - Rj\m)} ) , (23) 



where the same convention for the sums has been used 
as in Eq. ( pi] ) and we use (. . .) to denote Monte Carlo 
averages. The height of the structure factor evaluated 
at a reciprocal lattice vector scales with the size of the 
system. More specifically, we have 



5(Qi) =iVexp 



8^ 

3^0 



(u 2 ) 



(24) 



for the first Bragg peak in the triangular lattice, where 
(u 2 ) denotes the mean-squared fluctuations. Expression 
( pi| ) has the usual Debye- Waller form and requires the 
lines to fluctuate harmonically around their average po- 
sitions. 

We stated above that there should be no true long- 
range order in a system of finite thickness (/?). Conse- 
quently, there should not exist real Bragg peaks. In order 
to understand this point better, we compute (u 2 ) using 
elastic theory for a system of finite thickness. The main 
contribution comes from transverse fluctuations, given by 
the expression 



1 



(K(R) 



E 



T[l — cos(K • R)] 

CssK 2 + C44&z 



(25) 



Here, the two-dimensional K integral is over the Brillouin 
zone, with Kbz = V^r/ao, and the discrete /c z -vectors 
are given by k z = 2irn/L z . We ignore the dispersion 
of the tilt modulus for the moment. Separating off the 
k z = term and integrating over the remaining ones, we 
obtain 



-( U± (R)- MJ _(0)] 2 ) 



TK 



BZ 



T 



■ HR/a ) (26) 



47iyc 66 C44 2irc 6e L : 

in the limit R 3> do- The first term is the result for 
the transverse fluctuations in a bulk superconductor and 
contains only k z ^ modes. The second term is due 
to the k z — mode and thus corresponds to motion of 
the average position of a vortex line. In a sample of finite 
thickness there is always only quasi long-range order since 
the second term in Eq. ( |2q ) grows logarithmically with 
the distance R. This term, however, is only relevant if 
R > ao exp (c660qL z /T) , which grows exponentially with 
the thickness of the sample. In the limit L z — > oo, the 
second term of Eq. (|2^) vanishes and we recover the well 
known result of true long-range order in the Abrikosov 
lattice. 

We shall use the bosonic superfluid density as a mea- 
sure of the entanglement in the system. The super- 
fluid density is defined by the response of the system 
to a transverse gauge field, or, equivalently, to the mo- 
tion of the walls of the container. We therefore have to 
study a system described by the Hamiltonian H v , ob- 
tained from the original Hamiltonian by the transforma- 
tion p — > p — mv. The superfluid density is then given 
by 



Ps_ = J 

p mN dvdv 

t b d 2 



mN dvdv 



lnZ v (/3,A), 



(27) 



where F v and Z v represent the free energy and the par- 
tition function for the system with ppving walls. It can 
be shown that the latter is given byE3 



(28) 



Z v (/3,A) = Z (/3,A)exp — W-v 



where the winding vector W is defined as 



(29) 



This definition deserves an explanation. Since the upper 
ends (r = (3) are just a permutation of the lower ends 
(t = 0), the only paths which give a contribution to the 
winding number are those which leave the system on one 
side and return on another via the periodic boundary 
conditions. For these paths we have 



dr- 



~~dT 



R I (/3)-R J (0)+in, 



(30) 



where L is the linear size of the system and n S Z D is 
an integer vector. We are only considering D — 2 in this 
work, but the winding number definition of the superfluid 
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density is valid in any dimension D. Combining Eqs. 
and (|28|) we immediately arrive at 



(W 2 



p DA 2 PN' 



(31) 



which is an exact relation. It turns out, however, that it is 
numerically difficult to get good statistics for the winding 
number in large systems since it can only be changed by 
global moves involving on the order of \/~N world lines. 
For large systems we then made use of a simpler measure 
of quantum effects, which consists of simply counting the 
number of world lines which do not end on themselves: 



N e = number of lines with R;(/3) ^ Ri(0). 



(32) 



This quantity does not need any global moves to be 
changed and is therefore simple to measure for arbitrar- 
ily large systems. It is interesting to consider the relation 
between N e and p s . Obviously N e > is a necessary con- 
dition for p s > 0. There is no guarantee, however, that 
the condition is also sufficient. In order to arrive at a 
more quantitative relation between the two quantities, 
we can reason as follows: The quantities p s /p and N e /N 
are intensive and should not depend on the total number 
of particles in the system. If we assume that they are 
proportional to each other, it follows from Eq. ( |3l| ) that 



(W 2 ) oc iV e 



(33) 



which is reasonable, given that the vector W is the sum 
of N e randomly directed vectors. Unfortunately, we can- 
not offer a derivation of the proportionality constant in 
Eq. (j3jj) and thus have to rely on empirical results: We 
define the quantity p e according to 



p 



q(A, 13) N e 
DA 2 (3 N 



(34) 



and compute a (A, (3) by requiring p e — p s for small sys- 
tems where p s can be computed using the standard wind- 
ing number definition (|3l| ) . It turns out that for the range 
of parameters considered here, we can ignore the A and 
(3 dependence of a (A, (3) and use 



a(A,/3) =4.04. 



(35) 



We emphasize, however, that a(A,/3) does depend on 
both [3 and A in general. Fi g. B shows the superfluidity, 
measured according to Eqs. ( plf ) and (34), as a function 
of A for two different system sizes. The good agreement 
between p e and p s , particularly the sharp onset of p e at 
melting (see Fig. m, is a consequence of the strong first- 
order character of the transition, and makes us confident 
that p e provides a good estimate for the superfluid den- 
sity in the larger Bose systems. Nonetheless, we stress 
that p e is not a real order parameter since, in contrast 
to p s , it does not vanish in the lattice phase when taking 
the thermodynamic limit. 



The energy of the system, both for bosons and vortices, 
can be computed using the thermodynamic definition 



E = T 2 — lnZ(A,/3). 



(36) 



Due to the different temperature dependence of the pa- 
rameters A and (3 in the two cases, the resulting expres- 
sions for the energy are different. It is useful to define 
the two expectation values 



* = (££ 



(R. 



"D \2 

, m ) 



2Ah 



s 2 = (J2J2 tK ^ r ^)/x 



(37a) 



(37b) 



corresponding to the two terms of the action (|22|). The 
important point is the temperature dependence of the 
term Si. For bosons, we have r oc 1/Tg and. Si oc T B - 
The expression for the energy then becomesEij 



E B =T B {MN -Si)+T B S 2 , 



(38) 



where the two terms correspond to kinetic and potential 
energies, respectively. The quantum-mechanical kinetic 
energy is decreased by the fluctuations of the world lines, 
since the fluctuations correspond to a smearing of the 
wave function. 

In the case of vortices, the denominator of Si is given 
by 



T 

A 2 t oc — . 



(39) 



If we want to treat J 7 as a, model Hamiltonian, we should 
ignore its internal temperature dependence through the 
penetration depth A(T). In this case, £; and £o are con- 
stants and both the terms <Si and £2 are proportional to 
l/T. Computing the temperature derivative in Eq. ( |3^ ) 
we obtain 



(F) =T(Si-MN)+TS 2 , 



(40) 



where T is again the Ginzburg-Landau free energy func- 
tional. The expectation value of J 7 is a very useful 
quantity to measure, but does not represent the energy 
of the vortex system. As was pointed out by Hu and 
MacDonaldE-3, the internal temperature dependence of J 7 , 
through A(T), accounts for microscopic degrees of free- 
dom that give a significant contribution to the energy. 
We will show below how Eq. ( f40| ) can be combined with 
scaling properties of the London free energy to obtain the 
correct expression for the energy. 

Another advantage of doing the simulation at finite 
temperature is that it allows us to compute the excita- 
tion spectrum of the system. The idea is to analyze the 
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dynamic structure factor S(Q,iu) n ), which can be ob- 
tained as the partial Fourier transform of the retarded 
density correlator 



5(Q,r) 



1 

N 



(p(Q,r)p(-Q,0)) 



(41) 



The normalization is chosen so that 5(Q, r = 0) is the 
static structure factor of Eq. (^3|) . The standard method 
for obtaining the excitation spectrum from iS(Q,mj„) is 
to perform an analytic continuation to real frequencies 
using the maximum entropy method (see Ref. 64 for an 
overview). Here we use a simpler approach, similar in 
spirit to analytical work by Nelson and Seungc3. If we 
assume the structure factor at momentum Q to be de- 
termined by one quasi-particle excitation, we obtain the 
single-mode approximation 



C c 



(W r , 



r 



(42) 



QJ 



Here sq and Tq are the excitation energy and inverse 
lifetime, respectively, of the quasi-particles, and Cq is 
related to the static structure factor of the system. This 
form of the structure factor leads to an exponential de- 
cay of correlations along the world lines according to 
S(Q,t) oc exp(-eq|r|). 

In order to understand the physical implications of the 
single-mode approximation we consider the superfluid 
phase, where analytical results are available. Eq. ( E^ ) 
is consistent with the two standard sum rules, i.e., the 
fluctuation-dissipation theorem and the /-sum ruleE3, 
provided that 



Q 2 



2mS(Q) : 



C Q = 2e Q S(Q), 



(43a) 
(43b) 



where S(Q) is the static structure factor and all quan- 
tities depend only on the magnitude of the wave vec- 
tor Q since the supe rflui d is isotropic. The excitation 
spectrum eq in Eq. (43a) is the Bijl-Feynman approx- 
imation, which, according to the above derivation, is a 
direct consequence of assuming one quasi-particle to sat- 
isfy the standard sum rules. However, it is well known 
that the Bijl-Feynman approximation significantly over- 
estimates the energy of the excitations in the superfluid, 
and one excitation therefore does not exhaust the sum 
rules. Nonetheless, we can still assume the single-mode 
approximation to give a valid description of the dynamic 
structure factor for small frequencies, and use the form 
Eq. to fit the numerical data, as illustrated in Fig. [|. 
The inset shows the corresponding excitation spectrum, 
the Bijl-Feynman spectrum computed from the static 
structure factor, and the theoretical Bogoliubov result 



V(Q). 



(44) 



The agreement is very good for small Q but the Bijl- 
Feynman approximation seriously overestimates the ro- 
ton energy. We have not performed an extensive analysis 
of the systematic errors in our method and the results 
should therefore be considered as numerical estimates 
rather than the final result on the excitation spectra in 
2D Bose systems. 



B. Isobaric Monte Carlo 

A disadvantage of performing a simulation of a com- 
pressible system at constant density is that a first order 
transition will not be sharp. Rather, we will observe the 
coexistence of the two phases at the transition, compli- 
cating the interpretation of the data. As was shown in 
Sec. ||, vortices in an applied magnetic field correspond 
to bosons at fixed chemical potential, and a simulation 
using the grand canonical ensemble would solve the prob- 
lem with coexistence. It is, however, technically difficult 
to vary the number of world lines in the system. Further- 
more, by varying the number of world lines, we can only 
change the density in discrete steps, which is especially 
problematic for small systems. 

Since it is the density rather than the total number of 
particles which is important, we can instead choose to 
simulate the Bose system at constant applied pressure. 
The relevant thermodynamic potential is then the Gibbs 
free energy, which is defined according to 

Z(T, P, N) = e -G{T,PM)/T = J dAZ{T ^ A N) e ~PA/ Tj 

(45) 

in direct analogy to Eq. (|J). The three thermodynamic 
potentials F(T,A,N), il(T,A,fi), and G{T,P,N) all of- 
fer perfectly valid descriptions of the 2D Bose system, 
and the choice is therefore a matter of convenience. Since 
G(T,P,N) does not correspond directly to the experi- 
mental situation of a superconductor in an applied mag- 
netic field, we have to make sure that the jump in den- 
sity is independent on whether we use G(T, P, N) or 
f2(T, A,p). We therefore rewrite the free energy accord- 
ing to F(T,A,N) = Af(T,p), where p = N/A. The 
pressure and the chemical potential are then fixed by the 
conditions 



P 



/' : 



dj_ 

dp 
d£ 
dp 



f, 



(46a) 
(46b) 



We now consider two phases with free energy densities 
/i(T, pi) and f2(T,p2). The coexistence condition at 
fixed p (fi continuous) is given by 



A/ = /i - h = ^Ap, 



(47) 







where Ap = p\ — p%. The analogous condition at fixed 
applied pressure (G continuous) is given by 



A/ 



Ap 
Pi 



(P + h) 



Ap 

Pi 



(P + / 2 ) 



(48) 



Comparison with Eq. 
have 



|) shows that in both cases we 



A/ = ^A,, 
op 



(49) 



which is nothing but the Maxwell construction. The 
jump in density will therefore be the same irrespective 
of whether we carry out the simulation at constant P 
or constant p. In the language of vortices, this means 
that the magnetization jumps observed in the simulation 
are directly comparable to experiments. It is, however, 
not possible to deduce the strength of the applied mag- 
netic field from the pressure P, just as it is not possible 
to compute p from P. According to Eqs. (^) we have 
pp = P + /, which requires the knowledge of the free 
energy of the system. 

The implementation of the isobaric quantum Monte 
Carlo algorithm is similar in spirit to the classical case 
(see, e. g., Ref. |6f]). Some details, which are mostly of 
technical nature, are given in Appendix |^. 



V. NUMERICAL RESULTS FOR THE 2DBCL 

One of the advantages of using the Bose model for 
studying the vortex system numerically is that the al- 
gorithm can be tested against known results for bosons. 
We begin by studying the transition from a Wigner crys- 
tal to a superfluid in the 2D Bose Coulomb Liquid, i.e., 
with A = oo. Fig. |] shows the structure factor and the 
superfluid density, computed according to Eq. (|34|), as a 
function of A. We have used (3 = 300 and M = 100. The 
lattice disappears at A m « 0.062, in perfect agreement 
with the results by Ceperley and Magro using a different 
numerical technique for the same systemcJ. The transi- 
tion is very sharp: The two structure factors shown as 
insets correspond to A = 0.06205 and A = 0.06245, so 
that the relative change in A is less than 1%. Simulta- 
neously, the superfluidity rises from p e = to p e — p, 
showing that the system is either superfluid or a crystal. 
Since our method of measuring the superfluid density 
via p e is not well established, we want to make sure that 
we obtain the same results when measuring p s using the 
standard winding number definition. Fig. ^ shows the 
same transition as Fig. || for two smaller systems, where 
the winding number has been measured directly. The 
dotted vertical line indicates the position of the transi- 
tion for the larger systems. The curves for the structure 
factors cross exactly at this line, showing that the posi- 
tion of the transition does not depend on the size of the 
system. The curves for the superfluid densities cross very 



close to this line, making us confident that there is only 
one transition in the system. 

The transition from a Wigner crystal to a superfluid is 
a quantum phase transition, which is driven by a change 
in h rather than Tb- Since none of the two phases 
changes dramatically with temperature, the properties of 
the transition will also not depend strongly on temper- 
ature. This can be understood intuitively in the world 
line picture: The world lines have to be long enough, but 
not strictly infinitely long. There exists an elegant way 
to check whether the temperature used in the simulation 
is low enough to capture the ground state behavior: The 
free energy of the systems has to be continuous across 
the transition, which implies 



AF 



AE; 



pot 



T B AS = 0. 



(50) 



In particular, for the ground state we have AEkin — 
—AEp t since Tb = 0. Fig. [| shows the kinetic and po- 
tential energies for a system with AT — 81. The transition 
from a lattice to a superfluid is indicated by a sharp drop 
in the kinetic energy, about 7%, due to the fact that the 
particles become delocalized. This decrease in kinetic 
energy is compensated by an increase in the potential 
energy. The almost perfect match between kinetic and 
potential energies in Fig. ^ also shows that the systematic 
errors in the simulation due to the discretization of the 
imaginary time are comparable to the statistical error. 
Thus, we expect our results to hold for the limit r — > 0, 
corresponding to continuous lines, and for /? — > oo, cor- 
responding to zero Bose temperature. 

In Fig. we show the energy of the 2DBCL as a func- 
tion of the de Boer parameter A. The transition cor- 
responds to a small change in slope of this curve. The 
change is more clearly seen in the inset, where we have 
subtracted the slope of the crystalline phase. Strictly 
speaking, since we are simulating at a finite temperature, 
there is a small downward jump in energy in addition to 
the change of slope. The statistical errors are too large 
to allow for a reliable estimate of this jump, however. 



VI. NUMERICAL RESULTS ON VORTEX 
LATTICE MELTING 

The above results on the transition in the 2DBCL can 
be directly translated into results for vortex lattice melt- 
ing. Thus, Fig. ^ shows a first-order melting transition 
from a lattice into an entangled vortex liquid. As long as 
we use the approximation A = oo, the entire melting line 
is given by the value of A m , and we find 



B m (T) 



$0 4A^ £;£ A 2 

A 2 V3 T 2 



(51) 



This result is not applicable to the low field regime with 
ao > A, since the value of A m does depend on the range 
of the interaction in this regime. For superconductors 
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that are not too anisotropic, the Lindemann criterion 
produces the melting lineu, 



B m (T) 



^2 47FC i T 2 - 



(52) 



which is identical to Eq. (|5lJ) with ei — s So ln(ao/2y / 7r^) 
and A m cx <? L , We therefore expect the Bose model 
to give a good description of vortex lattice melting in 
YBCO, where the anisotropy, e « 1/7, is moderate. For 
more anisotropic superconductors, such as BSCCO, the 
electromagnetic coupling between pancake vortices be- 
comes important, an. [effect which is not accounted for 
by the Bose model£3't3. It is important to remember 
that the shape of the melting line is mainly determined 
by the temperature dependence of the characteristic line 
energy £o rather than by details of the simulation. We 
will therefore focus on the properties of the melting tran- 
sition itself, rather than on the shape of the melting line, 
in this work. 

The structure factor shown in Fig. || is related to the 
structure factor measured-in small-angle neutron scatter- 
ing (SANS) experiments^. The main difference is that 
our structure factor measures the correlation in the po- 
sition of the vortex cores, whereas the neutrons are scat- 
tered by the magnetic field of the vortices. This limits 
the resolution of SANS to length scales larger than the 
magnetic penetration depth A. 

The elastic theory for the Bose model of the v ortex 
lattice, i.e., with the elastic moduli given in Sec. II B, 



predicts (u 2 ) oc A. Using Eq. ( p4| ) it is possible to mea- 
sure the mean-squared fluctuations numerically, and the 
result is shown in Fig. ||. The dashed line is a fit to the 
data for small A and we find 



0.9615 Aal 



(53) 



As is seen from the fit, elastic theory works well for A < 
0.05. Closer to the melting transition the renormalization 
of the elastic moduli becomes important, and the linear 
dependence of the fluctuations on A breaks down. Just 
below the melting transition we find (u 2 ) « 0.065, leading 
to a Lindemann number of 



c L « 0.25. 



(54) 



This result is in agreement with usual estimates for vor- 
tex lattice meltingu. 

The amount of entanglement is always very small in the 
lattice: We typically find N e /N w 0.026 just below the 
transition, compared to N e /N « 0.56 in the liquid just 
above the transition. Entanglement therefore increases 
by a factor 20 at the transition already for a system with 
N = 81, and we have direct melting into an entangled 
vortex liquid. This is in .afffefiment with recent numerical 
results on either modelsc3'E2l, and with flux transformer 
experimental. 

In order to further characterize the melting transition 
we consider the energy, which was defined in Eq. (36). 



For the case of vortices, where Z — exp (— T/T), we 
obtain 



(55) 



The second term in this expression accounts for the inter- 
nal temperature dependence in the effective free energy 
and was neglected in early numerical work on vortex lat- 
tice melting, which lead to much too small estimates for 
the latent heat of the transition. 

Although it would be possible to measure the quantity 
in Eq. ( pq ) directly in a simulation, we adopt a simpler 
approach hcrdiS. By direct inspection of the free energy 
(0) we see that we can make use of the scaling form 



F{{Ri{r)},T) « £ {T)a F{{Ri{T)}), 



(56) 



where T is dimensionless and temperature independent. 
The approximation only neglects the temperature depen- 
dence of A(T) in the interaction, and Eq. (|56| ) therefore 
becomes exact in the limit A — * oo. Combining Eqs. (|55|) 
and @) we obtain the result 



E = {T)[l- 



T de 
7~o~dT 



i + t 2 



t 2 ' 



(57) 



where we have assumed £o{T) = £q(0)(1 — t 2 ) for the 
last equation. In the simulations, we always measure the 
quantity {T) and express all results in terms of 



/ 



NL z e (T) ' 



(58) 



where Tq is the GL free energy for a perfect lattice of 
lines. Thus, / measures the additional GL free energy, 
due to thermal fluctuations, per vortex and length in 
units of £q(T). The real energy of the system can then 
be easily computed using Eq. (|57|). In Fig. || we show the 
GL free energy of the vortex system as a function of A. 
The transition is clearly seen as a discontinuity in / with 
Af f» 0.013. Using Eq. (p7|), we compute the entropy 
per vortex and layer and obtain 



AS = 0.013 (1 + f 



£o(0)d 

) rp ) 



(59) 



where d is the layer separation. With parameters for 
YBCO (d w 12A, A Qb w 1400A, and T m « T c w 92AT) 
we obtain AS f=i|0.35fc_B, in rough agreement with exper- 
imental resultscJ. 

It is interesting to consider how much of the latent heat 
comes from the Josephson energy in our simulation, i.e., 
from the first term in Eq. . It has been shown numer- 
ically that approximately half of the latent heat comes 
from the Josephson energy in a model for a layered super- 
conductor, a relation which is expected to become exact 
for a line modeled. Here we reach the same conclusion in 
a slightly different way: Comparing Eqs. (pi) and (|40|) 
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we see that the bosonic kinetic energy is the negative of 
the Josephson energy. Since it was shown above that the 
changes in kinetic and potential energies for the quantum 
system have to cancel exactly, it follows that the change 
in the Josephson energy equals the change in the in-plane 
interaction energy. 

There exist a number of methods for testing numeri- 
cally whether a transition is first order. A widely used 
method is to make a histogram of the data for the en- 
ergy and look-^or two peaks, indicating the coexistence 
of two phasesEa. The method did not work very well in 
the present case, since the transition is too sharp: Co- 
existence of phases was only observed for small systems, 
where it is difficult to resolve the two peaks in the en- 
ergy distribution. The sharpness of the transition for 
large system is sufficient evidence that the transition is 
really first order. 

We now turn to the case of finite A, where the system 
has a finite compressibility. According to the Clausius- 
Clapeyron relation, 



As 



AB dH m 
'ItT dT 



(60) 



we now expect a change in the density of the vortex sys- 
tem at the transition. When going from the lattice to the 
liquid phase we have dH m /dT < and As > 0, leading 
to AB > 0. The vortex system therefore shows anoma- 
lous melting, just as ice, with the melted phase being 
the denser one. In the case of the vortex system, the 
reason for the anomalous behavior can be understood 
in an intuitive way: When the lines begin to entangle, 
the line tension will produce an attractive force between 
them, leading to an increased densitjO. If we assume 
H m w B m , which is true as long as A > ao, we can 
rewrite Eq. ( |60|) as 



A/ = 87rApA(T) 2 



(61) 



Fig. |lfj shows the jump in density for two different values 
of A. The interaction between the lines is smaller than 
for A = oo and the melting transition is shifted towards 
lower values of A m . The results are summarized in Ta- 
ble |. There is no dramatic change in the properties of 
the transition when going from A = oo to A rs ao. This 
shows that the logarithmic vortex interaction is a good 
approximation for most of the phase diagram. The agree- 
ment with the Clausius-Clapeyron relation in Eq. (pll) is 
very good for A = 1.06ao, proving that the approximation 
H m w B m is valid and that the transition is really first 
order. For A = 0.72a , the right hand side of Eq. (|6l] ) 
is too small, since \dB/dT\ < \dH/dT\. Using the fact 
that A/ does not depend strongly on A/ao, we can use 
the value A/ rs 0.013 and obtain a very simple expres- 
sion for the density change of the vortex lattice at the 
melting transition, 



A P k 5.2 x IQ-^/XiTf 



l.lxl0 6 [GA 2 ] , , 
= \(T) 2 (62) 



This expression is again in good a 
tization measurements on YBCOc 



jreement with magne- 



VII. PROPERTIES OF THE LINE LIQUID 

We now turn to the properties of the vortex liquid. 
It has been shown above that the vortex lattice melts 
directly into an entangled vortex liquid, which is equiv- 
alent to a bosonic superfluid. One important quantity 
in the vortex liquid is the correlation length l z in the 
direction parallel to the lines. If this correlation length 
grows very large close to the melting line, an apparent 
disentangled liquid could be observed in a sample of finite 
thickness. On the other hand, we want the correlation 
length to be larger than the layer distance in order to 
separate melting of a line system from decoupling. A lot 
of experimental effort has been invested into measuring 
l z , mainly using flux transformer techniques. The latest 
results indicate that there exists a shoxt but finite 2-axis 
correlation length in the vortex liquioQ. 

One of the nice features of the Bose model is that the 
question of the z-axis correlation length in the vortex 
system is related to the fundamental question of the roton 
minimum in a Bose system. In Sec. || we have shown that 
the coordinate z along the vortex lines corresponds to the 
imaginary time r for bosons. This relation can be made 
dimensionally correct by writing 



9 2 r 



2e l z 
T 



(63) 



According to Eq. [42] we expect the retarded density corre- 
lator to decay exponentially, with thexharacteristic time 
tq = H/eq. Using Eq. (|6^) we obtained 



UQ) 



2eo eq ' 



(64) 



The largest correlation length is found for the Q that 
minimizes eq, i.e., for the roton minimum. Since the 
roton minimum is located at Q ~ 27r/ao, this length is 
related to the single line correlation length. 

Fig. [ll] shows the excitation spectra for compressible 
and incompressible 2D Bose systems. At small Q, the 
compressible system has a phonon branch, whereas the 
excitations in the incompressible systems are plasmons. 
The roton minimum, however, does not change with the 
range of interaction. Since it is the roton minimum which 
determines the critical velocity of the superfluid, this sup- 
ports our claim that the 2DBCL does not differ qualita- 
tively from other 2D Bose liquids with regard to super- 
fluidity. The roton minimum collapses when the system 
solidifies, indicating 1) the appearance of a periodic lat- 
tice and 2) the disappearance of superfluidity. Note that 
the roton energy is large compared to the Bose temper- 
ature used in the simulation, Tg ~ 0.003 g 2 compared 
to Eroton ~ 0.026 g 2 . This explains why we observe the 
ground state behavior in the simulation. 
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Recalculating the roton minimum in terms of Eq. ( |63| ) 
we find 

This result is in good agreement with our expectations 
from elastic theory and scaling theory: The relevant 
length scale in the lattice is ao, which should be rescaled 
by anisotropy. The length l z is much larger than the 
layer spacing for YBCO, proving that the transition is 
really due to melting of continuous lines. As soon as the 
lattice is formed, l z becomes infinite and entanglement 
disappears. 
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VIII. CONCLUSIONS 

We have used the mapping to 2D bosons for studying 
the vortex system in a high-T c superconductor numeri- 
cally. An advantage of this approach is that it allows us 
to perform the simulation without introducing an artifi- 
cial lattice, which would cut off small scale fluctuations 
and introduce spurious pinning effects. Furthermore, by 
comparing the numerical results to known properties of 
the corresponding Bose system, we are able to control the 
systematic errors inherent in any numerical simulation. 

We find a first-order vortex lattice melting transition 
into an entangled vortex liquid. The results for the la- 
tent heat and the magnetization jump agree well with 
experiments for moderately anisotropic superconductors 
such as YBCO. We do not expect such good agreement 
for strongly anisotropic materials, since the Bose model 
describes vortices in terms of elastic lines, rather than 
pancake vortices, and ignores the electromagnetic cou- 
pling between the layers. 

In the language of bosons, we have found a sharp tran- 
sition from a Wigner crystal to a superfluid even for the 
case of a logarithmic interaction (the 2DBCL). Since the 
latter system does not have a Bose-Einstein condensate 
even in the ground state, this shows that the relation 
between superfluidity and Bosc-Einstcin condensation is 
rather subtle. It is amusing to note that the main result 
in this work can be predicted using two simple physical 
ideas: 1) the Landau quasi-particle picture and 2) the 
Bijl-Feynman form for the excitation spectrum in a Bose 
system. 

The emphasis of this work has been on the generic 
phase diagram for a system of interacting elastic lines. 
In this picture, the transition described above can be un- 
derstood in terms of a competition between "energy" and 
"entropy" : The ordered lattice always has the lowest "en- 
ergy" , but there are many more entangled states. The 
entanglement also offers a natural explanation for the 
anomalous melting characteristic of vortex lattice. In or- 
der for this scenario to work, however, the lines have to be 
long enough. Thus, one of the main results of this work 
is to elucidate the importance of the third dimension for 
vortex lattice melting. 



APPENDIX A: COMPUTING THE FOURIER 
SUM FOR THE INTERACTION 

We consider a rhombically shaped system with side L 
and the smaller angle 9. In the simulation we have used 
9 = 60°, but nothing is lost by solving the more general 
problem here. Since we want to be able to choose the 
range of the interaction A large compared to the size of 
the system, it is important to account for the periodic 
boundary conditions in the interaction. In the special 
case of the vortex interaction this is most easily achieved 
by solving the London equation using the correct periodic 
boundary conditions. We introduce a basis ei,e2 with 
e\-e 2 = cos 9. The system is then spanned by the vectors 
a, = Le^, arid the corresponding reciprocal lattice vectors 
are 



2tt 



(e 4 



(Al) 



for = (1,2), (2,1). The solution to the London equa- 
tion 



(1-A 2 A)G(R,A) = A 2 <S(R) 



is then given by 
G(R, A) 



2ttA 2 
L 2 sin 



^ 1 

G 



A 2 G 2 ' 



(A2) 



(A3) 



where the sum is over all reciprocal lattice vectors G = 
nibx + n 2 h 2 . This solution is proportional to the mag- 
netic field from one vortex line and in the case of L A, 
the function G(R, A) reduces to K (i?/A). Since there 
is no such simplification for gen eral L, we will have to 
evaluate the Fourier sum in Eq. (A3) numerically. It is, 
however, very time consuming to sum over n\ and n 2 di- 
rectly, and we have therefore decided to evaluate one of 
the sums analytically. This can be done with some effort 
using contour integration, and the resulting expression is 



G(R, A) 



2(1 + cos 9) 



^cos(At 1 )G JV (A,i 2 ), (A4) 



N 
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with t± t 2 = tt(R\ ± R2)/L and 



Cat(A, t 2 ) 



cosh[(7r/2- M)ajv(A)] 
aAr(A) sinh[7raAr(A)/2] 

sinh[(7r/2- |t 2 |)oAf(A)] 
oac(A) cosh[7raAr(A)/2] 



N 



even, 



N odd, 



where 



1 



2tt 2 A 2 1 + cos(9 1 



:N 2 



COS ( 



(A5) 



(A6) 



Eq. (A4) replaces an algebraically converging two- 



dimensional sum with an exponentially converging one- 
dimensional one, making the numerical evaluation trivial. 



APPENDIX B: SAMPLING PERMUTATIONS 

In order to capture the effects of Bose statistics, we 
have to allow the world lines to switch their endpoints. 
This can only be done by cutting out sufficiently long 
segments of a number of lines and trying different ways 
of connecting the loose ends. Updating the lines therefore 
becomes a two-step process: First we decide on how to 
connect the lower ends of the cut with the upper ones 
and then the new line segments are built up using the 
bisection method proposed by CeperleyE3. 

The problem of connecting the upper and lower ends 
of a cut of N world lines is difficult for the simple reason 
that there exist Nl ways of doing it. One way around the 
problem would be to choose only a small number of neigh- 
boring lines for the Monte Carlo updates, thus keeping 
N small. On the other hand, we know that superfluid- 
ity is due to large scale particle exchanges; the number 
of particles participating in a Monte Carlo update has 
to be large enough to change the winding number of the 
system. Different approaches have been suggested for 
solving this problem and some of them can be found in 
Ref. || 

The method we use for connecting the ends of the cut is 
a random walk in the space of permutations. The statis- 
tical weight of interconnecting lines is only determined 
by the kinetic part of the action and we introduce the 
weight 



W; 



exp ■ 



2A 2 A 



(Bl) 



for connecting the lower end of line i with the upper 
end of line j over a distance A in the imaginary time 
direction. The total weight of a permutation a of N 
lines is then given by 



JV 



W(a) = J] 



UK 



(B2) 



The probability of the same permutation is the weight 
normalized by the sum of the weights of all possible per- 
mutations, and this normalization factor is practically 
impossible to compute if N is large. The random walk 
algorithm selects a permutation in the following way: We 
first select how to connect the first line with a probabil- 
ity proportional to the weights. Thus, the probability of 
connecting the lower end of the first line with the upper 
end of line j is 



Pij 



Efc=l w lk 



(B3) 



After this we select how to connect the second line to any 
of the remaining N — 1 lines, again with the probability 
proportional to the weight. The procedure is repeated 
until all lines have been connected. The acceptance prob- 
ability for this permutation is then given by 



= n 



Ei=fc W kv(l) 



k=l 



l^l=k 



(B4) 



If the permutation is not accepted, the random walk is 
started again. The algorithm is exact as long as the 
identity permutation is the cheapest permutation, i.e., 
as long as A(a) < 1 for all possible permutations, some- 
thing which is almost always true even in the superfluid 
phase. It is still not possible to sample the permutations 
of arbitrarily large systems since the acceptance rate de- 
creases with N and it becomes necessary to try very many 
random walks before an acceptable permutation is found. 
However, since the weights lUy can be computed once and 
stored in a table, each random walk is very fast. Using 
this algorithm we were able to sample the permutations 
of N < 49 particle almost exactly. For larger systems, 
we typically restricted ourselves to N w 25 neighboring 
lines. 



APPENDIX C: ISOBARIC MONTE CARLO 
ALGORITHM 

In order to implement the isobaric Monte Carlo 
method, we first rewrite the discretized action ( p2|) using 
the system size L as a length scale, 

5[{xJ, L] — L 2 (X ^"' + 2 1 AV X "" )2 - MN ln(L 2 ) 

i,m 

+ £ tKq (^^i) +(3PL 2. (CI) 



V A 



Here we have added a term corresponding to the external 
pressure (cf. Eq. j45|)) and a logarithmic term coming 
from the change of variable, R = Lx, in the integral 



»=i 



M N 



M N 



nn^ 



3 MiV ln(L 2 ) 



/ II 1I' /J -'-— • ( C2 ) 
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Using Eq. (CI) we can compute the change in all ther- 
modynamic quantities as a function of the area of the 
system. The only problem is with the interaction term, 
which does not scale in a simple way with the area. In the 
present simulation, the interaction is stored in a lookup 
table and it would be too time consuming to reinitial- 
ize these tables every time the volume changes. We are 
helped by the fact that the volume changes are small, 
allowing us to use a Taylor expansion, 



U [L 2 {1 +e)] pb U(L 2 ) + L< 



dU 



1 



d 2 U 



-T 4 

2 rf(L 2 ) 2 ' 



where 



u(l 2 ) = tK o (HmA) ■ 



(C3) 



(C4) 



i<j,m 



The interaction and its two derivatives with respect to 
the area are evaluated for fixed L during the simulation. 
As long as the changes in L are small, we can us e th e 
Taylor expansion for evaluating the change in Eq. (CI). 



Ra the r than sampling the fluctuations in L, we used 
Eq. (CI) to optimize the action with respect to L af- 
ter the line configuration has been updated once. This 
gives the correct expectation value for the density and 
the other thermodynamic quantities. 
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superfluid 



normal fluid 



FIG. 1. Schematic phase diagram of the 2DBCL in the 
parameters A and f3. The solid lines indicate real phase 
transitions and quantum effects are important in the shaded 
region. For vortices in high-T c superconductors we have 
A 2 = T 2 /2eie al and (3 = 2e L z /T. The dashed line is 
obtained by keeping B (do) fixed and varying T. This line 
passes through all three phases for thin samples and is shifted 
to larger f3 as L z is increased. 
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FIG. 2. The superfluid density of the 2DBCL measured us- 
ing the winding number (dots) and the entanglement density 
(solid line) as explained in the text. We have used the param- 
eters /3 = 300, M = 100, and N = 25, 36. Both curves are 
fitted using the value a(A, 0) — 4.04, which was established 
by requiring the last point (A ~ 0.071) to have p e = p- 




FIG. 3. Fit of the single mode structure factor S(Q,iuj„) 
to measured data for a compressible system with N = 64, 
M = 100, and f3 = 300 at Q/Qi = 0.5. The fit is very good for 
low frequencies and fails for higher ones as expected. The in- 
set shows the excitation spectrum obtained from our method 
compared to the Bijl-Feynman and Bogolioubov spectra. 
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FIG. 4. Structure factor and superfluid density at the tran- 
sition from a Wigner crystal to a superfluid for the 2DBCL. 
We have used parameters f3 = 300, M = 100, and N = 64, 81, 
and have measured the superfluid density using p e as de- 
scribed in the text. The peaks in the structure factor dis- 
appear in a sharp transition at A ~ 0.062. Simultaneously, 
the superfluid density rises from p e ~ to p e ~ p. The tran- 
sition is very sharp, with a relative error A A/ A < 1%. In 
the language of vortices, the data shows a first-order vortex 
lattice melting transition into an entangled vortex liquid. 
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FIG. 5. The same transition as in Fig. ^ for smaller systems 
with N = 25, 36. The superfluid density is measured directly 
through the winding number. The vertical dotted line shows 
the position of the transition for the larger systems in Fig. ^[ 
The curves for the structure factor cross at this line, showing 
that the position of the transition is not shifted with the size 
of the system. The curves for the superfluid density also cross 
very close to the vertical line, although the statistical errors 
are larger. We have removed the error bars from some data 
points for clarity. 
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FIG. 6. Potential and kinetic energies at the transition 
from lattice to superfluid for the 2DBCL with N = 81 and 
/3 = 300. The mismatch in energy is less than 10% of the 
total jump, which is comparable to the statistical error. The 
relative change in kinetic energy at the transition is 7%. The 
potential energy has been shifted an arbitrary amount for 
clarity. 
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FIG. 7. Energy of the 2DBCL as a function of A for a 
system with N = 81, M = 100, and (3 = 300. The transition, 
indicated by the vertical dotted line, corresponds to a change 
in slope. This can be seen more clearly in the inset, where we 
plot E — cA and have chosen c to match the slope below the 
transition. 
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FIG. 8. Mean squared fluctuations of in the vortex system. 
The data points were obtained from the structure factor in 
Fig^. The dashed line is a fit to the data for small A and gives 
(u 2 ) ~ 0.9615 Aa 2 . The simple linear behavior breaks down 
for A > 0.05 due to renormalization of the elastic moduli. The 
Lindemann number computed from the size of the fluctuations 
just before melting is cl ~ 0.25. 
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FIG. 10. The change in density for two systems simulated 
with the isobaric Monte Carlo method. We have used the 
parameters /3 = 300, M = 100, and N = 64. The verti- 
cal dotted line indicates the position of the transition in the 
2DBCL. The transition is shifted to smaller values of A as a 
result of the weaker interaction. 
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FIG. 9. The GL free energy per particle and length for 
vortices with A = oo. The energy of a perfect vortex lattice 
has been subtracted. The melting transition is clearly seen as 
a discontinuity in the energy, Af ~ 0.013. The same systems 
as in Fig. ]4 were used. 
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FIG. 11. Comparison of the excitation spectra for com- 
pressible and incompressible systems. The parameters are 
13 = 300, M — 100, N = 64, and the compressible system has 
A ~ 1.06ao. The compressible system has a phonon branch 
for small Q. As the range of interaction is increased, the 
phonons turn into plasmons. The roton minimum is indepen- 
dent of the range of interaction, but collapses as the system 
crystallizes. 




TABLE I. Jump in energy and density for three systems 
with N = 64, /3 = 300, and M = 100 at different pressures. 
The relative error is less than 5% in the measured quantities 
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A/ao 


A m 


A/ 


87rA 2 Ap 


oo 


oo 


0.0622 


0.0127 




5 


1.06 


0.0596 


0.0113 


0.0116 


1 


0.72 


0.0567 


0.0091 


0.0085 
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